Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann 
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On-site boundary conditions are often desired for lattice Boltzmann simulations of fluid flow in 
complex geometries such as porous media or microfluidic devices. The possibility to specify the 
exact position of the boundary, independent of other simulation parameters, simplifies the analysis 
of the system. For practical applications it should allow to freely specify the direction of the flux, 
and it should be straight forward to implement in three dimensions. Furthermore, especially for 
parallelized solvers it is of great advantage if the boundary condition can be applied locally, involving 
only information available on the current lattice site. We meet this need by describing in detail how 
to transfer the approach suggested by Zou and He jj to a D3Q19 lattice. The boundary condition 
acts locally, is independent of the details of the relaxation process during collision and contains no 
artificial slip. In particular, the case of an on-site no-slip boundary condition is naturally included. 
We test the boundary condition in several setups and confirm that it is capable to accurately model 
the velocity field up to second order and does not contain any numerical slip. 

PACS numbers: 02.70.-c - Computational techniques; simulations 

47.11.-j - Computational methods in fluid dynamics 



I. INTRODUCTION 



The lattice Boltzmann method (LBM) is a widely used 
method for the simulation of fluid flow 2] . It solves the 
Boltzmann equation on a discrete lattice and it has been 
proven that the Navier Stokes equations can be recov- 
ered [1, H| ■ The method has been successfully applied to 
the simulation of flow in porous media 0] , colloidal sus- 
pensions [7H10j| , liqu id-g as phase transitions and multi- 
component flows [i"lT - ll3l ] spinodal decomposition [13, [Hj], 
and many more applications. 

In spite of the wide range of applications of the LBM 
there is sitll little consensus on how to implement bound- 
ary conditions in the LBM. For some applications, espe- 
cially for complex geometries in technical applications, 
rather simple approaches like on-site bounce back [l6[ 
rules are preferred [l7| , and on the other hand quite com- 
plex methods to implement exact boundary conditions 
have been proposed [l||, [l9| . A promising approach of 
velocity boundary conditions by Zou and He [l| for 2D 
simulations has been generalized to 3D with the restric- 
tion to the inflow being perpendicular to the boundary 
plane by Kutay et al. [5]. However, to our knowledge, 
a generalization to 3D with variable inflow direction has 
not yet been presented. Apart from this restriction, in 
the terms used in Ref. [5| some of the prefactors have 
to be revised (the correct ones can be found in Ref. [20j|). 
but for the application studied by Kutay et al. , the terms 
used in their work might be appropriate. However, es- 
pecially if the influx direction is not aligned with the 
computational lattice, our slightly more general expres- 
sions have to be used. We follow the ideas of Zou and 
He [l[ and derive flux boundary conditions with variable 
influx direction for a D3Q19-lattice [2l[ meaning that in 



three dimensions the velocity space contains 19 discrete 
vectors. Zou and He [l| have demonstrated a derivation 
for the D2Q9 model and shortly sketched the applica- 
tion to pressure boundaries in a D3Q15i model, where 
i stands for an incompressible model of equilibrium dis- 
tribution functions [22[ . Already Zou and He point out 
that besides the basic idea of applying a bounce back 
rule to the non-equilibrium part, a further modification 
is necessary to achieve the correct transversal momen- 
tum. A suitable choice for this correction depends on 
the lattice type. Zou and He give an expression for the 
D3Q15i lattice for the case of pressure boundaries. In the 
subsequent publications on D3Q19 lattices [H,[2(|, which 
is one of the most commonly used lattice types nowa- 
days, it is assumed that the flux direction is restricted 
to the direction normal to the boundary plane and that 
this symmetry is also reflected in the distribution func- 
tions on the boundary nodes. In our generalization we 
drop this restriction and consistently derive the transver- 
sal momentum corrections. We investigate the accuracy 
of this boundary condition and highlight the special case 
of on-site no-slip boundary conditions included in this 
approach by simply setting the velocity equal to zero. 
Examples for possible applications are microfluidic de- 
vices [23[, i.e., microscopic channel structures which are 
specially designed to modify a given flow profile by the 
roughness or wettability of the walls or the geometry of 
the channels. One example for those structures are mi- 
cromixers |24| . In general, if simulating such devices, it is 



not always possible to align all walls with the Cartesian 
planes. In these cases one needs a boundary condition 
which is capable to specify the velocity in an arbitrary 
direction, depending on the orientation of the channel to 
be simulated. 

One might also think of applications in porous media Q , 
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where flow through discretized samples of stones are sim- 
ulated. On the boundaries of the microscopic pores, a 
no-slip condition has to be applied. This is a special case 
of velocity boundary conditions with the velocity being 
zero. Since the channels cannot be aligned with the com- 
putational lattice, the question raises, how large the error 
introduced by the discretization is. The on-site velocity 
boundary conditions brought forward in the present pa- 
per can be used as a replacement of the bounce-back rule 
for the no-slip condition if the velocity is set to zero. In 
contrast to the usual bounce-back rule the position of the 
wall is independent of the BGK relaxation time. This fact 
is of great advantage when analyzing the permeability of 
a discretized sample. 

Although the assumption of the fluid velocity being 
zero on the boundaries does not hold for several cases 
in microfluidics, the no-slip condition is highly impor- 
tant for many cases. Therefore, a considerable effort 
of research has been spent to develop no-slip boundary 
conditions pi [l9l [25l |26| . Some approaches turned out to 
contain an artificial slip length depending on various de- 
tails of the simulation method, whereas other attempts 
involve non-local calculations like the evaluation of a ve- 
locity gradient to extrapolate the flow field beyond the 
boundary. In contrast to that, the boundary condition 
we propose is of local type and allows to specify the ve- 
locity on the node exactly with vanishing slip length. 
Further, the boundary condition is of great benefit for 
hybrid simulations, i.e., simulations in which two sim- 
ulation methods are coupled to simulate fluid flow [27i — 
[29j . The main goal of such hybrid simulations is to save 
computing time. A computationally cheap method is ap- 
plied to simulate the flow on a more coarse-grained level, 
whereas in some regions, where for example interactions 
on the atomistic level arc relevant, a different simulation 
method comprising more details is applied. The two sim- 
ulation methods are coupled for example by exchange of 
mass, momentum and energy between the two domains 
via their respective boundary conditions. One can think 
of different setups: an LB simulation can be embedded 
into a finite element based Navier Stokes solver and re- 
solve one region in more detail. Another example case 
is that in an LB simulation one region is resolved even 
on the molecular level by means of a Molecular Dynamics 
simulation. Practically, in current hybrid simulations the 
coupling is implemented within an overlapping region [30| 
of the two simulations, but it would be a great advance if 
one could manage simply to couple two boundary condi- 
tions without any overlap being needed. The possibility 
to generally determine the velocity on a boundary node 
in an LB simulation is one step towards this goal. 
The remainder of this paper is structured as follows: in 



the following section we describe the simulation method 
in general and introduce our notation of the lattice vec- 
tors. Then, we shortly review different boundary con- 
ditions in the literature in Sec. lIIII After that, we de- 
rive and discuss the velocity boundary condition for the 
D3Q19 lattice in Sec. lIVI We separately discuss the spe- 
cial case of the no-slip condition in Sec.fV] Numerical 
results are presented and discussed in Sec. lVII and finally, 
we draw a conclusion in the closing section of the current 
paper. 



II. SIMULATION METHOD 

The lattice Boltzmann method is a numerical method to 
solve the Boltzmann equation Eq. (JXJ) on a discrete lat- 
tice. The Boltzmann equation describes the dynamics of 
a gas from a microscopic point of view: in a gas, particles, 
each with velocities Vj, collide with a certain probability 
and exchange momentum among each other. For ideal 
collisions total momentum and energy are conserved in 
the collisions. The Boltzmann equation expresses how 
the probability /(x, v, t) of finding a particle with veloc- 
ity v at a position x and at time t evolves with time: 

vV x / + F-V p / + ^ = 0(/), (1) 

where F denotes an external body force, V xp the gradi- 
ent in position and momentum space, and f2(/) denotes 
the collision-operator. Bhatnagar, Gross, and Krook [3l| 
proposed the so-called BGK dynamics, where the colli- 
sion operator Q is chosen as a relaxation with a charac- 
teristic time r to the equilibrium distribution f( eq \v, p). 

&U) = -\(f-f (eq) ) ■ (2) 

The equilibrium distribution function for athermal mod- 
els depends on the local density p(x, t) and the veloc- 
ity field v(x, t). The lattice Boltzmann method [32[ dis- 
cretizes the probability density / in space and time. 
The discrete Boltzmann equation, which is solved by the 
LBM can be rigorously derived from the Boltzmann equa- 
tion (33j . The discretization, and especially the analytic 
expression for the equilibrium distribution depends 
on the lattice type. We use a D3Q19-lattice which is a 
very popular lattice type for 3D LB-simulations. On each 
lattice site 19 values /i(x, i) are stored, each of them as- 
signed to a lattice vector c^. We use the notation that 
the vectors are the i th column vector of the matrix 
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with the lattice speed of sound c s = for the D3Q19 
lattice and the lattice weights 
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The pressure p = c^p turns out to be proportional to the 
density and the dynamic shear viscosity is given by @, HH 
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III. BOUNDARY CONDITIONS 



Figure 1: The geometry of the D3Q19 lattice with lattice 
vectors c; as defined in Eq. (J3J). 



The geometry is shown in Fig.[TJ [49]. The local density 
at a lattice point can be obtained by summing up all fi, 
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pom) = ^2fi( x ,t), 



(4) 



and the streaming velocity is given by 
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v(x,i) 



(5) 



We express all quantities in lattice units, i.e., time is mea- 
sured in units of update intervals and length is measured 
in units of the lattice constant. For practical applications 
a suitable mapping to physical units based on a dimen- 
sional analysis has to be applied. 

In the lattice Boltzmann method two steps are performed 
in an alternating way: 

1. The "streaming step": propagate each of the dis- 
tribution functions /,; to the next lattice site in the 
direction of its assigned lattice vector Ci . 

2. The "collision step": on each lattice site relax the 
probability functions fi towards the equilibrium 
value f^ eq \v,p). In BGK dynamics this is accord- 



ing to Eq. ©. 



The equilibrium value is obtained by discretizing 

the Boltzmann distribution. Several expressions of dif- 
ferent order have been proposed, where we use the popu- 
lar form involving terms in the velocity up to the second 



On the boundary nodes, the distribution function as- 
signed to vectors Ci pointing out of the lattice move out 
of the computational domain in the propagation step, 
and the ones assigned to the opposing vectors are un- 
determined because there are no nodes which the distri- 
butions could come from. Therefore, on the boundary 
nodes, special rules have to be applied. 
These boundary conditions can be chosen in various man- 
ners. Periodic boundaries are realized by propagating the 
fi leaving the computational domain on the one bound- 
ary to the boundary nodes located on the opposite side 
of the domain. Closed boundaries are commonly im- 
plemented by a so-called mid-grid bounce-back rule [2j , 
which means that the distributions fi pointing out of the 
domain are copied to /•,-, for which Cj = — c^, i.e., locally, 
on each lattice site, the undetermined values are filled 
with the ones which would stream out of the domain 
without collision on the boundary node. They enter one 
time step later into the simulation domain again (39j . 
However, for many questions in fluid dynamics it is re- 
quired to determine the pressure or the velocity field at 
the boundary. The first is known as Dirichlet boundary 
condition, and the latter as Neumann boundary condi- 
tion. In the Neumann case the flux on the boundary 
of the domain is fixed, whereas in the Dirichlet case the 
pressure is given as a boundary condition. 
Zou and He [l[ have proposed how to implement Dirich- 
let and Neumann boundary conditions on a D2Q9 lattice 
and shortly sketched how to apply it for a D3Q15i simula- 
tion. Kutay et al. Q have transferred this proposal to a 
D3Q19 lattice. However, their approach is derived under 
the assumption that the in- and outflow velocity is al- 
ways perpendicular to the boundary plane, and oriented 
along one of the main lattice directions (cj, i — 1 . . . 6). 
We generalize this to inflow with arbitrary direction in 
SecHVl 
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Often more elaborated boundary conditions are ap- 
plied. Chen and co-workers [37[ and Ginzbourg and 
d'Humieres [l9| suggested extrapolation of the fi on the 
first and second layer of the lattice to the nodes outside 
the domain. These extrapolated values can be thought 
of as the lattice populations propagating into the do- 
main and arriving on the boundary nodes in the next 
streaming step. Inamuro and co-workers have introduced 
a counter-slip to compensate a numerical slip which oc- 
curs when applying on-site bounce-back [25j . Skordos 
came up with an approach where additional differential 
equations are solved on the boundary nodes to calcu- 
late the unknown populations (l8| . Ansumali and Karlin 
have developed a LB no-slip boundary condition from ki- 
netic theory [Z(|, and, more recently, d'Orazio et al. 
and Tang et al. [35| came up with thermal boundary 
conditions which also involve an extrapolation scheme 
and bounce-back with counter-slip respectively. Ladd 
and Verberg have developed a boundary condition with 
a resolution of the position of the wall on a sub-grid 
level, which is especially required if suspended particles 
are modeled [j| 0, E3| • Schiller and Diinweg [Sj use a 
reduced set of distribution functions on the boundary 
nodes. For their reduced D3Q19 model they derive equi- 
librium distributions and propose a multi relaxation time 
dynamics and a special collision operator on the bound- 
ary. 

Latt et al. have compared and discussed several of these 
approaches in Ref. [45(. They also include the boundary 
condition by Zou and He [l| in their discussion. As in- 
dicated by Latt et al. a generalization of the boundary 
conditions proposed by Zou and He is still not provided. 
However, a general local boundary rule which can be ap- 
plied in a simple way on each node separately, would be 
desirable. 

We derive such a boundary condition in the following 
section. Our generalization of the velocity boundary con- 
dition proposed in Ref. [l[ only involves the distribution 
functions defined on the local boundary node and allows 
by very simple and computationally cheap steps to set 
the velocity on the node to a distinct vector. The de- 
sired value is obtained exactly and we cannot detect any 
artifacts like a numerical slip length or bends in the ve- 
locity profile. 



IV. GENERAL ON-SITE VELOCITY 
BOUNDARY CONDITION 

As mentioned in the previous section, we extend the 
boundary condition by Zou and He [l[ to a D3Q19 lat- 
tice. We derive the boundary condition for the bottom 
plane (z = 0) in detail and give the results for the other 
planes in the appendix. They can be derived following 
the same steps. 

The boundary conditions are derived by using the set 
of equations consisting of Eq. (Q} and the components of 



Eq.©: 

pv x = h + h + fs + h + ho 

-(/2 + /11 + /12 + /13 + /14), (9) 

PVy = /3 + fl + fll + fib + fid 

-(/4 + /8 + /12 + /17 + /18), (10) 
pv z = fs + h + /is + fu + fn 

-(/6 + /10 + /14 + /16 + /18). (11) 

Due to the continuity relation §f + V • (pv) = 0, we are 
free to specify only three of the four variables (p and the 
three components of v) on the boundary. If we fix the 
tangential velocity v x , v y on the bottom-layer of the lat- 
tice, and the density to a given value po, the z-component 
of the inflow velocity v z can be calculated from Eq. (fTTj) 
and Eq. flU, 

Vz = l - — [/1 + /2 + /3+/4 
PO 

+ h + fs + fll + fl2 + /l9 (12) 
+ 2(/ 6 + / 10 + / 14 + / 16 + / 18 )] , 

where the /, pointing out of the system appear with a 
prefactor of 2, and all in-plane components appear with 
weight 1. The components pointing into the system, 
fbi /<), /i3j /i5j an d /17, which are undetermined after the 
streaming step, do not appear at all. With Eq. (fT2|) Neu- 
mann (or pressure) boundary conditions can be applied 
by specifying po on the boundary and using Eq. Q12[) to 
calculate v z . If Eq. (fT2|) is written in the form 

p = r^— Ifi + h + h + U 

l-v z 

+ /7 + /8 + /11 + /12 + /19 (13) 
+ 2(/ 6 + f 10 + fu + he + /is)] , 

all three components of the velocity can be specified and 
Eq. (Ti"3|) is used to calculate the density p. This is the 
Dirichlet case, or flux-boundary condition. Again, the 
undetermined populations /s, /g, /13, /15, and fu do not 
enter the calculation. 

We have used two out of four equations (Eqns. (|9"j)- (fTTj) 
and Eq. (0J), but we still have to compute the five /j 
pointing into the computing domain. Following Zou and 
He [l| we assume that on the boundary the bounce-back 
condition is still valid for the non-equilibrium part /* of 
the single particle distribution fi. 

ft =fi~ ft q) ■ (14) 

The bounce-back condition in +z-direction (in normal 
direction to the boundary) would read as 

n =/ 5 -/ 5 (e9) =/6-/ 6 (e9) =/6, (15) 
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which leads by taking f^ eq ^ and /g 6 ^ from Eq. (Q~4 



to 



same plane. Our Ansatz thus reads as follows: 



h = h ~ w 6 p 

+ w 5 p 
2w 5 



1 + ^ 



1L 

1 



= /e H — ^ = /e + 7/w z ■ 

C s 6 

Here we make use of the fact that the distribution func- 
tions in Eq. (J6j) are approximated by taking only terms up 
to 2 nd order in v into account. However, this approxima- 
tion could be applied directly to Eq. (fT6|) as well. For the 
derivation of the boundary condition it is needed, other- 
wise the higher order terms would introduce anisotropic 
effects in the boundary rule. 

Generally, in the collision step (in Eq. ©) higher order 
terms may be taken into account for the bulk, but for the 
boundary conditions a qualitatively different approach, 
like a higher order extrapolation scheme, has to be con- 
sidered when aiming for higher order accuracy. 

For the D3Q19 lattice, however, we need two more equa- 
tions. To keep the symmetry of the problem, we assume 
bounce-back of the non-equilibrium part for all popula- 
tions fi. This results in four equations, 





h 


= /l4H 


P ( 
6 


+ v x ) - K , 


(21) 




/l3 


- fio H 


g V z 


- v x ) + N* , 


(22) 




fib 


= AsH 


~ —{v z 
gl 


+ v y )-N*, 


(23) 


(16) 


fir 


= AeH 


Pi 
6 y 


- v v ) + N z y . 


(24) 



h 




2wg 

- ~7rP\ Vz - 


I" Vx) , 


(17) 


fl3 


= /ioH 


2^13 / 
- 2 P(V Z 


- Vx) , 


(18) 


fib 


= /«H 


2wi 5 

- 2 P( Vz 


+ Vy) , 


(19) 


fir 


= /ieH 


2w 17 

- 2 P(Vz 


~Vy), 


(20) 



so that the system of equations is overdetermined. 

Therefore, following the Ansatz by Zou and He |l| for 
the pressure boundary condition on a D3Q15i lattice, we 
introduce two new variables N£ and N*, the transver- 
sal momentum corrections on the z-boundary for distri- 
butions propagating in x and y-direction, respectively. 
These terms turn out to vanish in equilibrium, but they 
are non-zero, if velocity gradients are present, e.g., when 
shear flow is imposed by the particular choice of the 
boundary conditions. It turns out that these expressions 
appear again in the stress tensor. They reflect the fact 
that by imposing a transversal velocity component on 
the boundary, also stress is imposed to the system. The 
transversal momentum corrections involve the popula- 
tions propagating in the boundary plane in the update 
rule of the boundary condition. We add the terms to 
the right hand side and assume that the same expression 
with opposite sign is needed for two of the vectors in the 



The system of equations (j2T])-((2l]), together with Eq. (jH) 
and Eq. ([T0|) is now a closed system. By Eq. (|9]) and 
Eq. (fT0"|) we specify the tangential components of the ve- 
locity v x and v y , which do not need to be equal to zero 
in our approach. Inserting Eqns. (|2T]) -([24 |) into Eq. (|9]) 
and Eq. (fl0|) , gives an exact solution for N* 
respectively: 



. and N*, 



K = J[/i + /7 + / 8 -(/a + /n + /i2)] 



-,PVx 



N: 



[fa + h + hi - {fi + h + /12)] 



-pv y 



(25) 



(26) 



These transversal momentum corrections can be inserted 
into Eqns. (f2T]) - ([24| again and together with Eq. ([T6]l we 
find explicit expressions for all unknown populations. 
Note that in Eq. (|25p and (|26[) it is required to sum over 
all in-plane contributions to the velocity in x- and y- 
direction and the weights are consistent with the lattice 
weights of the fi appearing in the above expressions. As 
expected, N* and N* vanish (up to 2 nd order which is 
our precision within this derivation), if we set all fi to 
their equilibrium value. The results for the other planes 
are given in the appendix. 

A general form for all boundary planes can be written 
down by introducing the normal vector on the boundary 
n, the tangential vectors tj = C; — (cj • n)n, and the 
notation /_j denoting the direction to which a population 
is bounced back c_ s ; = — c,. From the populations fi 
assigned to a direction pointing into the wall, the new 
populations /_j in opposite direction can be calculated 
as 



f—i fi 
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■Ci-v-£ti-v+-^/j (t< • Cj) (1 - |cj • n|) . 

(27) 

Due to the particular choice of Eqns. (|2"T ]) -([2l) ]) or Eq. ([2"T]) 

respectively, it is possible to specify the velocity to an ex- 
act value on the lattice site. The rules presented here are 
independent on the relaxation rate in the collision step, 
since all calculations involve only the known values of 
the fi and equillibrium functions. Relaxation is calcu- 
lated separately after all unknown fi are calculated and 
the macroscopic velocity and density are preserved dur- 
ing collision. There are no restrictions on the orientation 
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of the inflow direction. Furthermore, all calculations are 
local on each lattice site. Apart from using only terms 
of first and second order in v for the equilibrium dis- 
tributions fl eq ^ in Eq. ([5]) no approximations are made. 
Therefore, we have derived a way to implement explicit 
local on-site boundary conditions which model the fluid 
field up to second order in the velocity. 
A similar scheme as ours has been proposed by Halli- 
day et al. [46[ for a D2Q9 lattice. These authors con- 
struct the unknown distributions locally on each lattice 
site starting from a Chapman-Enskoog analysis. During 
their derivation they have to choose a set of variables 
they consider as free variables. This is similar to the ap- 
proach of introducing the transversal momentum correc- 
tions in order to be able to solve the system of equations. 
Halliday et al. find results for the unknown populations 
which involve the components of the strain rate tensor 
calculated from the known populations. From this point 
of view it might be possible to apply a scheme similar to 
the one proposed by Halliday et al. to a D3Q19 lattice. 
However, in three dimensions the systems of equations, 
in the generality of Ref. (46|, might become difficult to 
handle. 

Special care has to be taken when connecting the in- 
and outflux boundary conditions at the corners and edges 
of the simulation domain with other types of boundary 
conditions that are applied on other boundary planes. If 
no-slip boundary conditions are assumed on the x- and 
y-boundary, one has to take care that the influx velocity 
tends to zero at the edges. We discuss the special case of 
no-slip boundaries as a subset of velocity boundaries in 
the following section. 



V. ON-SITE NO-SLIP BOUNDARY CONDITION 

The on-site velocity boundary condition proposed in this 
paper includes an important special case: setting the ve- 
locity v = results in a no-slip-boundary for non-moving 
boundaries. Therefore, this boundary condition can also 
be used as a replacement of the mid-grid bounce back 
rule. However, even more generally, moving boundaries, 
e.g., moving shear plates, can be implemented by im- 
posing the wall velocity v on the boundary nodes. The 
position of the wall is exactly on the lattice nodes. This 
is in contrast to most no-slip boundaries proposed in the 
literature, where the wall position is assumed at half the 
distance between two nodes. However, in many of those 
approaches the exact position of the wall depends on the 
BGK relaxation time. This is not the case for our ap- 
proach. 

One of the pillars of the LBM is local mass conserva- 
tion 45j, which should be fulfilled not only in the bulk, 
but also on closed boundaries. However, some extrap- 
olation schemes may be less accurate in this point [44| , 



whereas for our on-site approach mass conservation is 
strictly fulfilled at the closed walls. 

The transversal momentum corrections similar to those 



given in Eq. (|25|) and Eq. (|26|) are given in the appendix 
for each coordinate plane, and both velocity components 
in each of those planes. They are corrections to the on- 
site bounce back rule. With these corrections taken into 
account, the velocity is exactly zero on the node. For 
edge nodes with both boundaries being implemented as 
described here, we suggest to first apply the bounce back 
rule for all fa pointing out of the computational domain 
and then to calculate the transversal momentum correc- 
tion. On an edge node only one tangential vector along 
the edge can be used for ensuring no-slip. 
Consider for example the edge between the xy-plane and 
the yz-plane, where the y-axis forms the edge. Contri- 
butions in the boundary planes known after bounce back 
are f u f 7> and / 8 in the xy-plane and f 5 , f 15 , and f 17 in 
the j/2-plane. However, to ensure no-slip one can define 



N xz 
V 



h 



(28) 



The correction to the distributions fi with i = 7,8,15, 
and 17 then is N* z a ■ U which has to be added to the 

takes into ac- 



distributions. The prefactor in Eq. 
count that the remaining slip velocity after bounce back 
is distributed among four populations obtained from the 
bounce back rule. Similar expressions can be written 
down for each edge. A general expression for the modi- 
fied bounce back rule is 
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f-i = fi~ T fj (tj ' 



) 1- 



1 - 
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(2) 



(29) 

where and n^ 2 ) denote the two normal vectors on the 
two boundary planes meeting at the edge under consid- 
eration. 

On the edges and corners, apart from the incoming pop- 
ulations, there are so-called "buried links" [26j], i.e., lat- 
tice vectors c$ for which the opposing lattice vector c_; 
points out of the domain, as well. The two lattice vec- 
tors c' 1 ' 2 ' = ± (n' 1 ) — n^ 2 )) make up the buried link on 
an edge node. In the following, lattice vectors with sub- 
script, Cj, denote distinct vectors as defined in Eq. ([3]), 
whereas vectors with superscript, c^, denote vectors 
which belong to the buried links, and which depend on 
the normal vectors on the individual boundary planes. 
The distribution functions assigned to the buried links 
have to be assigned separately. We choose them such 
that they contribute to the same density according to 
their lattice weight: 



c(l,2) 



18 , 

hT.fi (i 

i=l v 



n (D x n (2) 



(30) 



n 



(i) 



n 



(2) 



V2 



Similarly, the distribution on the resting node is chosen 
as / 19 = ^a/t 1 ' 2 ' = 12 /(b2). The weights are always 
determined by the number of fi which contribute to the 
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sum and their respective lattice weights Wi according to 
Eq. ||7J| . In Eq. (|3"01) six fi with lattice weight yg and ten 
/i with weight ^ contribute, which makes up an overall 
contribution of 1§. To reduce this to the desired lattice 

36 

weight, we have to divide by 22, and to obtain a value 
for the resting node, we multiply by 12, because of the 
twelve times larger lattice weight of the resting node. 
At the edges surrounding the in- and outlet planes, on the 
other hand, one needs either pressure or velocity bound- 
ary conditions, depending on the boundary type used for 
in- or outlet. For velocity boundaries one has to take 
care that the velocity profile decays to zero, so the no- 
slip boundary condition just described can be used on all 
edges. For pressure boundaries one prescribes a density 
p = po which we can be used to calculate the distribution 
assigned to the buried link: 



(1,2) 



22/(1.2) 



14 



(31) 



where fig is then calculated as fig = 12 f^ 1 '^ . 
According to Maier et al. [26| no-slip boundaries cannot 
be enforced on convex boundary nodes. However, slip 
along the edge can be reduced by correcting all distribu- 
tions fi traveling into the interior of the system by 



Ni 



1 
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(2) 



(32) 



which follows the same idea as Eq. (|29|) and 
Eqns. (f2"T|) - (|2"4"|) : momentum in a direction parallel 
to the surface, which would remain on a node after 
applying the boundary rule, is removed by modifying 
those populations that will afterwards propagate back 
into the bulk of the system. In principle, one could 
split Eq. ([29)) into two steps: first, apply bounce back 
for all populations leaving the system, and then correct 
the populations traveling away from the edge by the 
term given in Eq. (|3"21 . For convex edges these are the 
populations traveling into the bulk, and for concave 
edges they propagate in the boundary planes. This 
opens a possibility to implement all rules in a single 
procedure, for which the normal vector n is stored on 
each lattice site by an integer number. The vector is 
obtained from the matrix M defined in Eq. ^ . For 
values between 1 and 6 Eq. (|2"7) is applied, for values 
between 7 and 18 cither Eq. (1321) applies or depending 
on the values stored on the neighboring nodes, the 
bounce-back rules corrected according to Eq. (1291) may 
be applied instead. This information, which expresses 
if the edge is concave or convex, can be obtained once, 
when the lattice is generated and may be stored in 
the sign of the lattice vector index. The normal vector 
points into the bulk and indicates the direction of the 
symmetry plane on the edge nodes. 

Finally, on the corner nodes one can define three normal 

vectors on the boundary planes meeting there, 

and n^ 3 ) . Similar to the buried links, there is a complete 



plane in which six vectors are located, that only couple 

to the simulation in the collision step. The buried vectors 
c (i...6) are ftie oneg £ or which C W . ( n (l) _|_ n (2) _|_ n (3)~j _ 

. Since the normal vector on this plane is not con- 
tained in the set of vectors for the D3Q19 lattice addi- 
tional indices are needed to mark the corner nodes. After 
bouncing back the known fi, the distributions assigned 
to buried vectors are set to 



.6) 
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i( 2 ) 



,(3) 
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(33) 
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18 

/in is set to 



if velocity boundaries are applied, or to /C 1 — 6 ) 
/(1...6) pressure boundaries are chosen. / 19 
^2 j(i- -6) or 12 f( 1 — 6 ) 1 respectively. A correction similar 
to Eq. (f2"9"|) is not necessary on the corner nodes[50j. 
In complex geometries there are points in which edges 
(convex or concave ones) meet planes which are oriented 
perpendicular to the direction of the edge. There, we 
propose to use bounce back for those populations which 
would leave the computational domain and to assign 
an appropriate value to the resting node, similar as de- 
scribed for the corner nodes. There are no buried links, 
because those links are located inside the boundary plane 
which the edge connects to, so the resting node must be 
set to fig = || X)j=i fi ■ I n total there are 6 planes and 
4 possible orientations of the edges, each of them either 
convex or concave, making up 48 more cases. However, 
since there are no buried links and a momentum correc- 
tion is not necessary either, the rules can be implemented 
easily in only a few lines of code. In the following section 
we show the results of tests of the boundary condition in 
simple geometries like Poiseuille flow between two plates, 
where the exact solution is known. As an example for 
more complex geometries we simulate the flow through 
a rectangular channel, where also edges and corners are 
involved. 



VI. NUMERICAL RESULTS 

We test our boundary condition by simulating a 
Poiseuille flow through a tilted channel. The size of the 
computational domain is 64 x 8 x 128 LB nodes, where 
the channel has a width of 20 nodes and is tilted by an 
angle a = arctan(^) s» 17.48°. This angle is chosen 
such that both ends of the channel intersect the xy-plane 
at the top and the bottom of the computational domain 
and that there are two lattice sites of wall at the left and 
at the right of the channel at the bottom and the top 
plane respectively. The flow through our test channel is 
simulated in three dimensions. However, for convenience, 
the y-direction is periodic. The walls (only in this test) are 
implemented as simple bounce-back nodes. Here we ap- 
ply the boundary conditions derived in Sec. lIVI as in- and 
outflux conditions and compare them to other implemen- 
tations. 
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We choose this simple test because the analytical solu- 
tion for the flow field is known and so we can estimate the 
numerical error. Usually, one would avoid to have walls 
not aligned with the computational lattice because of the 
staircase like discretization of the walls, which brings an 
additional discretization error into the simulation. This 
discretization error can be avoided by simply aligning 
the channel with the computational lattice. However, if 
more complex structures, like, e.g., Y-channels for ap- 
plications in microfluidics are simulated, it may happen 
that always at least one channel is not aligned with one 
of the Cartesian directions. A technical workaround, if 
appropriate boundary conditions are missing, is to sim- 
ulate a very long channel so that in the first section of 
the channel, the fluid can relax to a steady flow profile, 
and only afterwards enters the actual simulation domain. 
However, this causes the computational effort to increase 
substantially. 

Knowing the width of the channel, the center of the in- 
and outlet, and a given velocity vo on the center line, one 
can calculate a Poiseuille flow field inside the channel, 



v P (x) 



vo 



x - x - 7(2: - Zp) 
Ax 



(34) 



where xq and zq denote the center of the simulation space, 
Aa; is the half width of the channel measured along the x- 
direction, and 7 denotes the increment due to the tilting 
angle, which is related to the components of the velocity 
by 7 = tana = £ = J« . 

We simulate flow through such a channel and apply dif- 
ferent in- and outlet boundary conditions. We use a re- 
laxation time t = 1 in all simulations presented here. 
However, we have checked that the results do not de- 
pend on this particular choice. After 5000 time steps a 
steady flow field is reached. However, to be sure that 
the simulations have converged, we simulate 20000 time 
steps until we evaluate data. 

To visualize the difference between simulation and theo- 
retical prediction we subtract the velocity on each lattice 
node and draw the resulting vector field in Fig.[3J The 
value of the velocity is scaled by a factor of 1500 for 
drawing the arrows. The colors are assigned the absolute 
value of the velocity after scaling the difference field. 
In Fig.[5]a) we apply the boundary condition used by Ku- 
tay et al. [5| to a case, where the restriction of the in- and 
outflow velocity parallel to the ^-direction introduces an 
error in the region close to the boundary. Note that in 
Ref. apart from assuming the velocity perpendicular 
to the boundary, the authors underestimate the transver- 
sal components, which may be of no importance in this 
case. We use the correct coefficients as presented very re- 
cently in Ref. [20j . but keep the restriction to the inflow 
perpendicular to the boundary, which has a much larger 
influence on the flow field. Not only the first and sec- 
ond layer of nodes close to the boundary are affected, 
but the boundary condition introduces vortices which 
have approximately the size of the diameter of the chan- 




Figurc 2: Velocity difference fields for different approaches: 
In- and outflow velocity constraint to the direction perpen- 
dicular to the boundary plane (a), In- and outflow velocity 
tilted and of parabolic shape as in the analytic solution, but 
with N£ and Ny set equal to zero (b), the boundary condi- 
tions derived in Sec. lIVI with the correct choice for N£ and 
Ny. difference field (c). The velocity difference vectors are 
scaled for drawing by a factor of 1500 and the absolute value 
of the velocity difference is color coded from blue (small) to 
red (large). 



nel. Therefore, the first step to generalize the bound- 
ary condition from Ref. Q to a case where the in- and 
outflow velocity has an arbitrary orientation, is to use 
Eqns. (fT71) -(|20 |) . as used in the simulation for which the 
result is shown in Fig.[3]b). It is obvious that this bound- 
ary condition still introduces vortices close to the in- and 
outflow. The strength of the vortices is smaller com- 
pared to the case shown in Fig.[2]a). However, the size 
of the vortices is comparable to the width of the channel 
here as well. The value of the tangential velocity on the 
boundary nodes differs from the value one inserts into the 
equations. By introducing the transversal momentum 
corrections N£ and N* in Eqns. (j2"Tj) - (l2"4"l) , the vortices 
disappear and the velocity takes exactly the value one 
specifies with Eqns. ([9l)-(fTTj) as one can see in Fig.[2]c). 
The remaining differencefield can be mostly ascribed to 
the discretization on the lattice. Each single step of the 
wall discretized to individual steps can be found in the 
flow profile. However, at the in- and outflux boundary no 
additional artifacts can be seen, which demonstrates the 
strength of our boundary condition. The velocity on the 
boundary nodes takes exactly the value which we specify, 
and therefore, no vortices are generated. 
The transversal momentum corrections iVJ and N* could 
also be understood in terms of a counter-slip similar to 
the approach of Inamuro et al. [25j ] , but the Ansatz how 
to obtain the unknown populations fi is different: we as- 
sume a bounce back rule for the non-equilibrium part 
of the distributions and end up with a linear correc- 
tion to the reflected populations, whereas the authors of 
Ref. 12 511 construct the unknown distributions based on 
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kinetic theory where the correction appears not on the 
level of the distribution functions but as a counter-slip 
on the level of the wall velocity. The values for the den- 
sity and the velocity inserted into the equilibrium dis- 
tributions in Inamuro's method are different from the 
ones used for bulk nodes. In our approach, however, the 
boundary nodes are treated similar as the bulk nodes: 
the velocity on the boundary node can be calculated by 
inserting Eq. ([25]) and Eq. ([36]) into Eqns. (j2T])-{l24j) and 
the obtained distributions together with the one from 
Eq. ([11]) and the density from Eq. (TJ3]) into Eq. ([5]). It 
turns out that the velocity calculated from Eq. ([5]) is ex- 
actly the one which is imposed at the boundary node by 
Eqns.(pT])-f24]). 

In all simulations we kept the tilting angle of the chan- 
nel constant, because the error of our boundary conditon 
is angle independent. We can quantify the quality of 
the boundary conditions by computing the ratio of the 
absolute value of the difference held and the calculated 
velocity field on each node. The obtained values are av- 
eraged over the first twenty layers of LB nodes from the 
boundary. 



| v(x) — v p (x) 
|v P (x)| 



dV, 



(35) 



Where the volume V contains those layers of lattice 
nodes, which are at most a distance of the channel width 
apart from the boundary of the computational domain. 
This captures approximately the vortices and provides a 
measure for the quality of the boundary condition. The 
results for the different cases shown in Fig. [2] are listed 
in the following tabular: 



Boundary condition 


relative error £ 


on-site velocity (Fig.[2]c) 


0.0996 


AfJ and N* set to zero (Fig.|2]b) 


0.126 


v x = v y = (Fig.[2]a) 


0.175 



Good agreement with the expected Poiseuille flow pro- 
file (Eq. ([M]) ) is reflected in small relative errors. Large 
numbers indicate deviations in the area, where the fluid 
fields are compared. We ascribe the remaining deviations 
to the discretization error of the wall and the accompany- 
ing uncertainty in the exact wall position in the present 
case of the staircase like discretization. We check this 
by increasing the resolution of the simulation by a factor 
of two. As we expect, the numerical error due to the 
staircase-like discretization decreases roughly by a factor 
of two to £ = 0.051. This shows that the staircase-like 
discretization introduces a first order error. Therefore, 
we need further investigations to see the second order 
accuracy of the in- and outflux boundary condition. 
As another test for the flux boundary condtion we sim- 
ulate a straight channel aligned parallel to the computa- 
tional lattice, again in a 64 x 8 x 128-domain with the 
same boundary conditions as for the inclined channel. 
The remaining relative error decreases to 0.00235, which 
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Figure 3: Poiseuille flow between two parallel walls driven 
by a body force. The simulation data are averaged in each 
lattice plane parallel to the walls and agree up to floating 
point precision with the calculated Poiseuille profile. 



is typical for Poiseuille flow simulations at this resolution 
in combination with a mid-grid bounce back rule on the 
boundary. 

We can further measure the quality of our boundary con- 
dition in a shear simulation. On a 32 3 lattice we apply 
periodic boundaries in x- and y-direction and impose a 
shear velocity of v x = ±0.02 with opposite sign on the 
top and bottom plane. We obtain a linear flow profile 
within floating point precision. There are no notable 
jumps between the first and second layer of LB nodes, 
which confirms that the strain rate tensor II is set up 
correctly on the boundary nodes. 

In a next step we simulate Poiseuille flow again, but this 
time we use a 32 3 lattice with periodic boundaries in 
y- and z- direction. We apply a body force [HI, 53| by 
adding a force term 



Av = — 

P 



(36) 



to the velocity in Eq. © in the whole simulation volume. 
The Poiseuille profile we expect is of the form 




X ~ X{) 



\,r 



(37) 



where the viscosity is given by Eq. ([5]). The velocity pro- 
file found in the simulation together with the expected 
Poiseuille profile is plotted in Fig. [3] The parabola con- 
tains no fit parameters. The velocity is exactly zero on 
the boundary nodes, whereas with a simple bounce-back 
a numerical slip can be observed, which results in a ve- 
locity of 4 x 10~ 5 for the same simulation setup without 
using the transversal momentum corrections N£ and N* . 
We have carried out this test with r = 1, but the data 
presented in Fig. [3] is obtained with r = 2 to ensure that 
our boundary conditions are not restricted to the special 
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case of t = 1. Apart from the influence of r on the vis- 
cosity (Eq. (|SJ)), our simulation results are not affected 
by the relaxation time. In particular, we do not see any 
r-dependent (numerical) slip. 

Also in this second-order test we find that the numerical 
error is of the size of the floating point precision on the 
computer. This underlines that our boundary condition 
reproduces the velocity field up to second order. 
Finally, we test our implementation by simulating flow 
through a square channel. The analytical solution for 
the velocity of a pressure driven flow in a b x b square 
channel is 48] 



v(x,y) 
C{x) 



~2r] 



C{x) 



, with 



r 3 

n=0 



(-iy 



cosh I - — — | cos 



(38) 

^ (2n+l)7T^ 



(2n + l) 3 cosh 



^ (2n+l)7T ^ 



where x E [—6/2,6/2] and y E [—6/2,6/2] are the coor- 
dinates in the cross-section, with the origin in the center 
of the channel. The pressure gradient Vp is imposed by 
pressure boundaries and the dynamic viscosity is known 
from Eq. (JU). The infinite sum in Eq. (l38l) can be trun- 
cated when a given accuracy is reached. We sum up 50 
terms and compare this approximation to our numeri- 
cal results on a 32 3 domain, i.e., b = 15 plus one layer 
of boundary nodes. In Fig.[4]a) we compare the analyt- 
ical solution from Eq. Q38[) with our simulation results. 
The velocity in z-direction is averaged over the y- and 
z-direction and the averaged value is plotted against the 
position in x-direction. A very good agreement of the nu- 
merical result with the analytical solution can be seen. 
For comparison, results for the node based bounce back 
rule are shown. For this boundary condition it is known 
that it is only first oderer accurate, which can be seen 
in the kink in the velocity profile close to the boundary 
nodes. It can however be made second order accurate 
by choosing the position of the wall somewhere (depend- 
ing on the BGK relaxation time) in between two nodes, 
which is known as the mid-grid bounce back Q • As one 
can see in Fig.@]a), if the wall-position is chosen correctly 
for the bounce-back rule, a satisfiying accuracy can be 
achieved, too (top-down triangles). Note that the posi- 
tion of the wall is shifted by half a lattice unit due to the 
different approach at the wall. 

In Fig.|4]b) the relative error depending on the size of 
the simulation is studied. There are three different errors 
involved. The truncation errors of the sum in Eq. (1381) can 
not be seen in this figure. If the sum is truncated after 
just a few terms, the error increases on two of the edges. 
In the corners of the simulation domain the errors due 
to the discretization on the lattice remains. This is what 
determines the accuracy for relatively small simulations. 
However, if the lattice is refined, this error decreases. It 
decreases with the square of the lattice constant which is 
typical for second order accurate schemes. Another error 
which dominates for large systems is the floating point 
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Figure 4: a) Velocity profile in a square channel aver- 
aged along the j/z-planes for the noslip boundary condition 
(squares), for the node based bounce back rule (triangles), 
mid-link bounce back (top down triangles), and the analyti- 
cal solution. The noslip boundary condition collapses with the 
analytical solution, whereas the on-site bounce back bound- 
ary condition shows a kink in the profile close to the boundary 
nodes. The numerical results are obtained on a 32 3 lattice and 
for the analytical solution the sum in Eq. (|38|l is truncated af- 
ter 50 terms. A 2D profile of the ^-velocity in a a;y-cross 
section is displayed as an inset. 

b) The relative error for system sizes between 8 and 64 nodes: 
in the corners next to the boundary the largest relative de- 
viations occur. Note that on the boundary nodes, where the 
truncation error of Eq. (I38|l is largest, the deviation between 
simulation and approximated analytical solution is negligible. 
Therefore, the deviation can be taken as a measure for the 
quality of the numerical result. 

c) The averaged error for different lattice resolutions versus 
the number of lattice nodes in each dimension confirming the 
second order accuracy of the no-slip boundary condition. 
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precision which is reflected in noisy data in the center of 
the simulation domain. This error is independent on the 
lattice constant. In Fig.|4]b) the relative error is shown for 
different system sizes. In the corners the error decreases 
with the system size, whereas the noise in the center is 
independent by the system size. In Fig.|4]c) we plot the 
mean error averaged over the whole system against the 
number of lattice nodes used for computation in each 
dimension. The slope of approximately 2 confirms the 
second order accuracy of the boundary condition. For 
the simulations presented in Fig. 2] pressure boundaries 
according to Eq. (IT21 are used and on the walls and edges 
we apply no-slip conditions as described in Sec.fV] For 
this plot we use only the range in which the lattice size 
dependent error dominates. For 64 lattice nodes in each 
dimension, the floating point precision in one of our post 
processing steps dominates the overall error. Therefore, 
we only use the smaller systems for this investigation. 
The inset in Fig. 2] shows the analytical velocity profile 
in a cross section perpendicular to the extension of the 
square channel. 



VII. CONCLUSION 

We have derived an explicit local on-site flux boundary 
condition for LB simulations on a D3Q19 lattice. Veloc- 
ity terms up to second order enter the derivation and this 
accuracy is also confirmed in the numerical tests. The 
in- and outflux velocity underlies no restrictions to any 
peculiar direction. We have demonstrated the numeri- 
cal accuracy by comparing simulation results for a flow 
through a tilted channel with the theoretical expectation 
of a Poiseuille flow. Remaining errors can be assigned 
to the discretization on the lattice and to rounding er- 
rors due to the floating point representation. We have 
tested the boundary condition in simulation of Poiseuille 
flow between two planar walls and in shear flow. In those 
tests the simulation data fits exactly the analytical solu- 
tions without any slip-parameter and independent on the 
BGK relaxation time. For this test we have used no-slip 
boundary conditions which are a special case included 
in the general velocity boundary conditions. Finally, we 
have tested the boundary condition by simulating the 
flow through a square channel. The scaling of the nu- 
merical error with the lattice resolution again confirms 
the second order accuracy. 
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Appendix 

Here we give the expressions for the other boundaries 
not treated explicitly in the text. We start with the top- 
plane where we implement outflux for our simulations. 
We obtain 

p = ~~t t/i + h + h + h 

v z + l 

+fy + fu + fvi + fa + fi9 + (39) 
2(/ 5 + h + /l3 + /is + /it)] 

or 

Vz = -1+-L/1+./2 + /3 + /4 
PO 

+fr + hi + /12 + fa + fw + (40) 

2(/ 5 + /9 + /l3 + /l5 + /l 7 )] 

with v z defined in positive z-direction. Here, the unde- 
termined populations after the streaming step are 
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with N§ and N* defined as previously, in Eq. ([25} and 
Eq. (|26j). 

For the left, right, front and back boundaries, which we 
do not use in this work one finds the following expres- 
sions. For the left (x — 0) boundary, 
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with 



and 



K = ^ [/3 + /15 + /16 - (.A + /17 + /is)] 
1 

—^PVy, 

K = ^[/ 5 + /l7 + /l5-(/6 + /l6 + /l 8 )] 

-\pv z . 

At the right boundary we have 
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At the front (y = 0) boundary, one hnds 
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At the back the density is given by 
1 



[A + k + k + k 



Vy + 1 

+ k + AO + A3 + Zl4 + /l9 

+ 2(/ 3 + k + fu + As + fie)} , 



or the velocity reads 



1 



= -1 + -LA + .A + A + /6 

Po 

+ k + flO + fl3 + fl4 + fl9 

+ 2(/ 3 + k + fu + As + Ae)] , 



and the distributions are 



(69) 



(70) 



(71) 



(72) 



k 


= k - T^PVy , 


(73) 


fl2 


= f 7 + ^(-Vy-V x ) + Ny, 


(74) 


k 


= fu + £(-v y + v x )-NV, 


(75) 


fl8 


= fu + ^(-v y -v z ) + Ny, 


(76) 


fl7 


= fw + ^(-v y + v z )-N*. 


(77) 
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